%SuffN24MONTHLY.m
% Model solution for Gold model

clear
tic

% Calibration
% Data targets for y10 and y1 are computed below based on annual averages
% Model moments here are set (manually) to match
% 1) meanr = mean of data 1 year rate first (will still change)
% 2) stdr, ar1r &sigm are moved manually for the model to match moments
%       for y10: std, ar1 and average (check mean y1 and adjust meanr)
% 3) gam is set to match average sensitivity 50/50 -- moved manually

meanr=   .000665;       
stdr =   .00349 ;       
ar1r =   .9877  ;

gam = 1.00079;   %Price trend growth:  max 1.011^(1/12)

bet = (1/1.035^(1/12));  

lam = 1/120;     % 1 over maturity of bond: 10 years, 120 months
vbar = 1;

N = 100;       %No of states

sigm = .0687 ;  %innov std of sdf
                %selected manually to match 10Y average yield 

%% Calibrating the model
% Get Markov process approximating a VAR: PI, rv 
%  rv  = [r ]    log rate
%

F = ar1r;
varrinno = (stdr^2)*(1-F^2);
F0 = 0;

% Tauchen program
bandwidth= 3.5;   
[PI,  rv] = tauchen(F, F0, varrinno, N, bandwidth);
    
% Spread component of M matrix - "log factors" -- see DerivationsSUFF.tex
y1v = meanr+rv;
xM = zeros(N,N);
b = -sigm;
for k=1:N
    eN = (y1v-PI(k,:)*y1v)';
    ek = (1/ (PI(k,:)*(eN.^2)'))^.5;
    eN = ek*eN;
    a = log(PI(k,:)*exp(-b*eN)');
    xM(k,:) = exp(-a - b*eN);
end

M = exp(-y1v).*ones(N,N).*xM;

%% Compute implied yields for 1 year and 10 years
ylam10Yv =   log (ones(N,1)./( ( eye(N) -(1-lam)*(PI.*M))\exp(-y1v) ) -lam +1 ) ;
lam1Y=1/12;
ylam1Yv = log (ones(N,1)./( ( eye(N) -(1-lam1Y)*(PI.*M))\exp(-y1v) ) -lam1Y +1 ) ;

% Unconditional probabilities row vector
eigvecu = null(PI'-eye(N));
piu = eigvecu'./sum(eigvecu)';

meany1M  =   piu*y1v;
meanylam1Y = piu*ylam1Yv;
meanylam10Y =   piu*ylam10Yv;

stdy1M =   (piu*((y1v).^2)-meany1M^2)^.5;
stdylam1Y = (piu*(ylam1Yv.^2)-meanylam1Y^2)^.5;
stdylam10Y =   (piu*(ylam10Yv.^2)-meanylam10Y^2)^.5;

ar1y1M =     ((piu*((y1v).*(PI*(y1v)))) - meany1M^2) /stdy1M^2 ;
ar1ylam1Y =  ((piu*(ylam1Yv.*(PI*ylam1Yv))) - meanylam1Y^2) /stdylam1Y^2 ;
ar1ylam10Y = ((piu*(ylam10Yv.*(PI*ylam10Yv))) - meanylam10Y^2) /stdylam10Y^2 ;

% Interest rate processes
calmod1M =  [meany1M    stdy1M     ar1y1M ]';
calmod1Y =  [meanylam1Y stdylam1Y ar1ylam1Y]';
calmod1Y =round(calmod1Y,5);
calmod10Y = [meanylam10Y stdylam10Y  ar1ylam10Y]';
calmod10Y =round(calmod10Y,5);
calmodmom = {'mean';'stddev';'AR1'};
Tabcalmod = table(calmodmom,calmod1Y, calmod10Y);

%% Solve for gold price

% Solved by Bellman iteration

vbarv = vbar*ones(N,1);
lastmaxdi=SuffNDiscCheckf(PI,M,vbarv,bet,gam);
if lastmaxdi>1
   return
end

[pv,pIv,vv,jstar,bfail]=SuffNSolBellf(PI,M,vbarv,bet,gam);
if bfail==1
        'Bellman equation iteration did not converge'
        return
end

pvRA= (eye(N)-((PI.*M)*(gam)))\vbarv;
    
futr1MRA = ones(N,1)-vbarv./pvRA;
fut1Mr = (((PI.*M)*gam)*pv)./pv;
fut12Mr = ((((PI.*M)*gam)^12)*pv)./pv;
  
pbar = vbar/(1-bet*gam);

sensv =   -log(pv(2:N)./pv(1:N-1))./(ylam10Yv(2:N)-ylam10Yv(1:N-1));
sensRAv = -log(pvRA(2:N)./pvRA(1:N-1))./(ylam10Yv(2:N)-ylam10Yv(1:N-1));

% Compact computation of variance of first difference

dPV = (repmat(log(pv),1,N)'-repmat(log(pv),1,N));

dPI = PI.*(repmat(piu,N,1))';
Ep2 = dPI(:)'*(dPV(:).^2);
Ep  = dPI(:)'*(dPV(:));
stdp = (Ep2-Ep^2)^.5;
% other variables added below
ModImpstat = [  piu*pv/pbar  stdp piu(jstar+1:N)*ones(N-jstar,1) ];
 
%% Plots Model Properties
figure(1)
subplot(2,2,1)
p=plot(12*ylam10Yv,pv/pbar);
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
grid
title('Gold prices vs 10-year rates')
subplot(2,2,2)
p=plot(log(pv/pbar),piu);
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
title('Unconditional probabilities of log gold price')
subplot(2,2,3)
p=plot( 12*ylam10Yv(1:N-1),sensv/12,'-');
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
grid
title('Sensitivity: -dln(Gold price)/dy10')
xlabel('10-year rate')
subplot(2,2,4)
p=plot(12*ylam10Yv,piu);
p(1).Color = [0, 0, 0]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
title('Unconditional probabilities of log 10-year rate')
xlabel('10-year rate')

%  exportgraphics(gcf,'figmod1.png','Resolution',300)


%% Model with exchange rate: Calibrate RER and Solve model for gold price 
% pst (pstar) gold price in foreign currency units

Nsti = 12;       % No of points for each dimension r & s
                 % For exploration run around 12-25, for final 100 (this is
                 % slow 30)
Nst = Nsti^2;

%% Calibrating the RER model and solve it
% Add exchange rate change std and correlation with interest rate innovation
% to Markov process
% Get Markov process approximating a VAR: PI, rv, sv
%  rsv  = [r s ]   - log rate and log rer change

% Empirical counterparts computed below, copied here manually
stds = .013;    % std of innv rer (log)
rhors = -.24;   % corr r_innov and s change (log)
                % Note: simulated correlations of s change with y1M and y1Y
                % are essentially identical. Can set rhors to corr with
                % y1Y.
                
Fst = zeros(2,2);
Fst0= zeros(2,1);

Fst(1,1)=ar1r;
varrinnost(1,1) = (stdr^2)*(1-F(1,1)^2);

varrinnost(2,2) = stds^2;
varrinnost(1,2) =rhors*(varrinnost(1,1)*varrinnost(2,2))^.5;
varrinnost(2,1)=varrinnost(1,2);

% Tauchen program 
bandwidth= 3.5;   
[PIst,  rsv] = tauchen(Fst, Fst0, varrinnost, Nsti, bandwidth);

% Spread component of M matrix - "log factors" 
y1vst = meanr+rsv(:,1);
xMst = zeros(Nst,Nst);
b = -sigm;
for k=1:Nst
    eNst = (y1vst-PIst(k,:)*y1vst)';
    ekst = (1/ (PIst(k,:)*(eNst.^2)'))^.5;
    eNst = ekst*eNst;
    a = log(PIst(k,:)*exp(-b*eNst)');
    xMst(k,:) = exp(-a - b*eNst);
end

Mst = exp(-y1vst).*ones(Nst,Nst).*xMst;

S = repmat(exp(rsv(:,2))',Nst,1);
% Set E(1+ds)=1 
sbarv= 1 - (PIst*exp(rsv(:,2)));
Sbarv = repmat(sbarv,1,Nst);
S = S+Sbarv;

%% Compute implied yields
ylam10Yvst = log (ones(Nst,1)./( ( eye(Nst) -(1-lam)*(PIst.*Mst))\exp(-y1vst) ) -lam +1 ) ;
lam1Yst=1/12;
ylam1Yvst = log (ones(Nst,1)./( ( eye(Nst) -(1-lam1Yst)*(PIst.*Mst))\exp(-y1vst) ) -lam1Yst +1 ) ;

tst = clock;  

vbarvst = vbar*ones(Nst,1);
[pstv,pIstv,vvst,bfailst]=SuffNSolBellstf(PIst,Mst,S,vbarvst,bet,gam);
if bfailst==1
        'Bellman equation iteration did not converge'
        return
end

timesolveFxmodel=etime(clock, tst);

%% Compare model to data

% GPIA10 Y10IA begm10 Y1IA10        from 79.10 to 23.12 all in LEVELs
% Furat1M10 Furat1Y10               from 79.10 to 23.12
% USDIndex in levels directly from FRED from 79.10 to 23.12
load datforSuff24
load datfutforSuff24MON

USDI  = xlsread('FRED-RTWEXBGS.xls', 1, 'F6:F536');

Tdat  = max(size(GPIA10));

tgam = gam.^[0:1:(Tdat-1)]';

pmod   =   interp1(ylam10Yv,pv,log(1+Y10IA/100)/12,'linear');
vmod   =   interp1(ylam10Yv,vv,log(1+Y10IA/100)/12,'linear');  

pmodst =   interp1(ylam10Yvst(1:Nsti,1),pstv(1:Nsti,1),log(1+Y10IA/100)/12,'linear');
pmodRA =   interp1(ylam10Yv,pvRA,log(1+Y10IA/100)/12,'linear');
fut1Mrmod10 = interp1(ylam10Yv,fut1Mr,log(1+Y10IA/100)/12,'linear');
fut1Mrmod1 =  interp1(ylam1Yv,fut1Mr,log(1+Y1IA10/100)/12,'linear');
fut1Mrmod10RA = interp1(ylam10Yv,futr1MRA,log(1+Y10IA/100)/12,'linear');
fut1Mrmod1RA =  interp1(ylam1Yv,futr1MRA,log(1+Y1IA10/100)/12,'linear');

fut12Mrmod10 = interp1(ylam10Yv,fut12Mr,log(1+Y10IA/100)/12,'linear');
fut12Mrmod1 =  interp1(ylam1Yv, fut12Mr,log(1+Y1IA10/100)/12,'linear');

figure(2)
p=plot(begm10(1:Tdat), GPIA10/mean(GPIA10),...
       begm10(1:Tdat), (tgam.*pmod/mean(tgam.*pmod)),...
       begm10(1:Tdat), (tgam.*pmodst./USDI)/mean((tgam.*pmodst./USDI)) )  ;
grid
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.55;
p(2).Color = [153 0 0]/255; 
p(2).LineWidth = 1.5;
p(3).Color = [1, 31, 91]/255; 
p(3).LineWidth = 1.5;
dateFormat = 12;
datetick('x',dateFormat)
title(['Real Gold Price: Model and Data ' ],'FontSize',14)
axis([datenum(begm10(1)) datenum(begm10(Tdat)) 0 2.5])
legend('Data', 'Model Yields only','Model Yields and USD','Location','NorthWest','FontSize',11) 

%   exportgraphics(gcf,'figmod2.png','Resolution',300)

figure(3)
p=plot(begm10(1:Tdat),[Furat1M10,fut1Mrmod1,fut1Mrmod10 ]);
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
p(2).Color = [153 0 0]/255; %The Penn red color code for the Penn Quakers logo is Pantone: PMS 201 C, Hex Color: #990000
p(2).LineWidth = .75;
p(3).Color = [153 0 0]/255; %The Penn red color code for the Penn Quakers logo is Pantone: PMS 201 C, Hex Color: #990000
p(3).LineWidth = 1.2;
title('Futures (1-month) to spot (with interest) ratio: Model vs Data')
legend('Futures ratio', '1-year real rate','10-year real rate' ,'Location','NorthEast') 
legend('Data', 'Model fits 1YRR','Model fits 10YRR','Location','NorthEast') 
dateFormat = 12;
datetick('x',dateFormat)
axis([begm10(1) begm10(Tdat)  .99 1.004])
grid

%  exportgraphics(gcf,'figmod3.png','Resolution',300)

figure(4)
p=plot(begm10(1:Tdat),[Furat1M10,fut1Mrmod1RA,fut1Mrmod10RA ]);
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
p(2).Color =  [128 128 128]/255;  % Gray / #808080
p(2).LineWidth = .75;
p(3).Color =[128 128 128]/255;  % Gray / #808080
p(3).LineWidth = 1;
title('Futures (1-month) to spot (with interest) ratio: Representative Agent Model vs Data')
legend('Futures ratio', '1-year real rate','10-year real rate' ,'Location','NorthEast') 
legend('Data', 'RA model fits 1YRR','RA model fits 10YRR' ,'Location','NorthEast') 
dateFormat = 12;
datetick('x',dateFormat)
axis([begm10(1) begm10(Tdat)  .99 1.004])
grid

%  exportgraphics(gcf,'figmod4.png','Resolution',300)

%% Calibration target empirical moments
% Calibrate gam by running annual growth rate regression 
% Set cutoff            50%   33%
cuto10 = .023;     %  .023;  .0153;
cuto1 =  .014;    % -.0024 .016  --- not used in paper

y10cal  = log(1+Y10IA(13:Tdat)/100);
dy10cal = log(1+Y10IA(13:Tdat)/100)-log(1+Y10IA(1:Tdat-12)/100);

y1cal  = log(1+Y1IA10(13:Tdat)/100);
dy1cal = log(1+Y1IA10(13:Tdat)/100)-log(1+Y1IA10(1:Tdat-12)/100);

samp10L = (y10cal<cuto10);
samp1L =  (y1cal<cuto1);

% Select with or without exchange rates
pmodfxcal = tgam.*pmodst.*(1./USDI) ;
pmodcal = tgam.*pmod ;
pmodRAcal = tgam.*pmodRA ;

dpdatcal = log(GPIA10(13:Tdat))  - log(GPIA10(1:Tdat-12));
dpmodcal = log(pmodcal(13:Tdat)) - log(pmodcal(1:Tdat-12));
dpmodcalRA = log(pmodRAcal(13:Tdat)) - log(pmodRAcal(1:Tdat-12));

% 10-year rate
stats = regstats2NW12(dpdatcal(samp10L),[dy10cal(samp10L) ],'linear');  beta1c=stats.beta; pval1c=stats.hac.pval; rsq1c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcal(samp10L),[dy10cal(samp10L) ],'linear');  beta2c=stats.beta; pval2c=stats.hac.pval; rsq2c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcalRA(samp10L),[dy10cal(samp10L) ],'linear'); beta3c=stats.beta; pval3c=stats.hac.pval; rsq3c=stats.adjrsquare; 

stats =   regstats2NW12(dpdatcal((1-samp10L)==1),[dy10cal((1-samp10L)==1) ],'linear'); beta4c=stats.beta; pval4c=stats.hac.pval; rsq4c=stats.adjrsquare; 
stats =   regstats2NW12(dpmodcal((1-samp10L)==1),[dy10cal((1-samp10L)==1) ],'linear'); beta5c=stats.beta; pval5c=stats.hac.pval; rsq5c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcalRA((1-samp10L)==1),[dy10cal((1-samp10L)==1) ],'linear'); beta6c=stats.beta; pval6c=stats.hac.pval; rsq6c=stats.adjrsquare; 

sensilow = [ beta1c(2) beta2c(2) beta3c(2)]';
sensihigh = [ beta4c(2) beta5c(2) beta6c(2)]';
sensilocuto = [  sum(samp10L)/(Tdat-11) nan nan]';
sensi = { 'Data' ; 'Model'; 'RA-model'};
Table5Sensi = table(sensi, sensilow, sensihigh, sensilocuto);

calgam ={'Sens low (model-data)';'Sens high (model-data)'; 'sum, Target=0'};
calgamval= [ (beta2c(2)-beta1c(2))  (beta5c(2)-beta4c(2)) (beta2c(2)-beta1c(2))+(beta5c(2)-beta4c(2)) ]';
Tabfitgam = table(calgam,calgamval);

% 1-year rate
stats = regstats2NW12(dpdatcal(samp1L),[dy1cal(samp1L) ],'linear');   beta7c=stats.beta; pval7c=stats.hac.pval; rsq7c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcal(samp1L),[dy1cal(samp1L) ],'linear');   beta8c=stats.beta; pval8c=stats.hac.pval; rsq8c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcalRA(samp1L),[dy1cal(samp1L) ],'linear'); beta9c=stats.beta; pval9c=stats.hac.pval; rsq9c=stats.adjrsquare; 

stats =   regstats2NW12(dpdatcal((1-samp1L)==1),[dy1cal((1-samp1L)==1) ],'linear'); beta10c=stats.beta; pval10c=stats.hac.pval; rsq10c=stats.adjrsquare; 
stats =   regstats2NW12(dpmodcal((1-samp1L)==1),[dy1cal((1-samp1L)==1) ],'linear'); beta11c=stats.beta; pval11c=stats.hac.pval; rsq11c=stats.adjrsquare; 
stats = regstats2NW12(dpmodcalRA((1-samp1L)==1),[dy1cal((1-samp1L)==1) ],'linear'); beta12c=stats.beta; pval12c=stats.hac.pval; rsq12c=stats.adjrsquare; 

% disp(['Regression: price change on 1-year yield, annual % change, DATA and MODEL'])
  [ beta7c(2) pval7c(2)  rsq7c sum(samp1L)/(Tdat-11) 
    beta8c(2) pval8c(2)  rsq8c sum(samp1L)/(Tdat-11) 
    beta9c(2) pval9c(2)  rsq9c sum(samp1L)/(Tdat-11)   ]  ;

  [ beta10c(2) pval10c(2)  rsq10c sum(1-samp1L)/(Tdat-11) 
    beta11c(2) pval11c(2)  rsq11c sum(1-samp1L)/(Tdat-11) 
    beta12c(2) pval12c(2)  rsq12c sum(1-samp1L)/(Tdat-11)   ]  ;

ModImp = {' E(p/pbar)' 'Std(dlnp)'  'Prob(I-out)' 'Sens low' 'Sens high' };
ModImpstat = [ ModImpstat beta2c(2) beta5c(2)  ];
Table4ModImp=table(ModImp);
Table4ModImpstat =table(ModImpstat);

% Interest rate calibration targets
% HERE everything is monthly

y10IAm = log(1+Y10IA/100)/12;
y1IAm =  log(1+Y1IA10/100)/12;
iUSDIm = -log(USDI);

stats = regstats2(y1IAm(2:Tdat), [y1IAm(1:Tdat-1)],'linear');  beta2=stats.beta; pval2=stats.hac.pval; rsq2=stats.adjrsquare; res2= stats.r;
stats = regstats2(y10IAm(2:Tdat),[y10IAm(1:Tdat-1)],'linear'); beta3=stats.beta; pval3=stats.hac.pval; rsq3=stats.adjrsquare; res3= stats.r;
stats = regstats2(iUSDIm(3:Tdat)-iUSDIm(2:Tdat-1),[iUSDIm(2:Tdat-1)-iUSDIm(1:Tdat-2)],'linear'); beta4=stats.beta; pval4=stats.hac.pval; rsq4=stats.adjrsquare; res4= stats.r;

cal1Y =  [mean(y1IAm)  std(y1IAm)   beta2(2) ]';
cal1Y = round(cal1Y,5);
cal10Y = [mean(y10IAm) std(y10IAm)  beta3(2)]';
cal10Y = round(cal10Y,5);
calmom = {'mean';'stddev';'AR1'};
Tabcaldat = table(calmom, cal1Y, cal10Y);
Table3bottomCalmom = table(Tabcaldat,Tabcalmod);

calpar    = {'mean y1';'innov stddev';'AR1';'risk price'};
calparval = [ meanr varrinno^.5 ar1r -sigm ]';
Table3topCalpar = table(calpar,calparval);

Table3topCalpar
Table3bottomCalmom
Tabfitgam

Table4ModImp
Table4ModImpstat

Table5Sensi

% Moments for exchange rate model (IID random walk): need stds and rhors (corr)
calirernams = {'std_innov_log_ds';'corr_innv_log_y1Y-ds'};
calirer1 = [std(iUSDIm(2:Tdat)-iUSDIm(1:Tdat-1)) corr(res2,iUSDIm(2:Tdat)-iUSDIm(1:Tdat-1))]';
TabRerTargets = table(calirernams, calirer1)


%% Analyze model deviations from data
% Measured in terms of log changes 12-month overlapping 

dpdat  =   log(GPIA10(1+12:end)./GPIA10(1:end-12));
dpmod  =   log(pmodcal(1+12:end)./pmodcal(1:end-12));
dpmodfx =  log(pmodfxcal(1+12:end)./pmodfxcal(1:end-12));
dpmodRA =  log(pmodRAcal(1+12:end)./pmodRAcal(1:end-12));

begfit1 = datenum(2004,1,1);    % 20 years 
endfit1 = datenum(2023,12,31);
begm1fit1 = begfit1<=begm10 & begm10<=endfit1;

begfit2 = datenum(2007,1,1);    % Best visual later period 
begm1fit2 = begfit2<=begm10 & begm10<=endfit1;

% Mean Squared Errors (normlized by total sum of square deviations from mean): NMSE
residdmod = dpdat-dpmod;
residdmodfx = dpdat-dpmodfx;
residdmodRA = dpdat-dpmodRA;

% 4 samples: full(f), late 20 years(lat), early(ear), late from 2007(lat2)
NMSEf(1,1) = sum(residdmod .^2)/sum((dpdat-mean(dpdat)).^2);
NMSEf(2,1) = sum(residdmodfx.^2)/sum((dpdat-mean(dpdat)).^2);
NMSEf(3,1) = sum(residdmodRA.^2)/sum((dpdat-mean(dpdat)).^2);
NMSEf=round(1-NMSEf,2);

NMSElat(1,1) = sum(residdmod(begm1fit1(13:end)).^2)  /sum((dpdat(begm1fit1(13:end))-mean(dpdat(begm1fit1(13:end)))).^2);
NMSElat(2,1) = sum(residdmodfx(begm1fit1(13:end)).^2)/sum((dpdat(begm1fit1(13:end))-mean(dpdat(begm1fit1(13:end)))).^2);
NMSElat(3,1) = sum(residdmodRA(begm1fit1(13:end)).^2)/sum((dpdat(begm1fit1(13:end))-mean(dpdat(begm1fit1(13:end)))).^2);
NMSElat=round(1-NMSElat,2);

NMSEear(1,1) = sum(residdmod(begm1fit1(13:end)==0).^2)  /sum((dpdat(begm1fit1(13:end)==0)-mean(dpdat(begm1fit1(13:end)==0))).^2);
NMSEear(2,1) = sum(residdmodfx(begm1fit1(13:end)==0).^2)/sum((dpdat(begm1fit1(13:end)==0)-mean(dpdat(begm1fit1(13:end)==0))).^2);
NMSEear(3,1) = sum(residdmodRA(begm1fit1(13:end)==0).^2)/sum((dpdat(begm1fit1(13:end)==0)-mean(dpdat(begm1fit1(13:end)==0))).^2);
NMSEear=round(1-NMSEear,2);

NMSElat2(1,1) = sum(residdmod(begm1fit2(13:end)).^2)/sum((dpdat(begm1fit2(13:end))-mean(dpdat(begm1fit2(13:end)))).^2);
NMSElat2(2,1) = sum(residdmodfx(begm1fit2(13:end)).^2)/sum((dpdat(begm1fit2(13:end))-mean(dpdat(begm1fit2(13:end)))).^2);
NMSElat2(3,1) = sum(residdmodRA(begm1fit2(13:end)).^2)/sum((dpdat(begm1fit2(13:end))-mean(dpdat(begm1fit2(13:end)))).^2);
NMSElat2=round(1-NMSElat2,2);

NMSEvar = {'yields only';  'yields &  rer'; 'RA yields'};
Table6topNMSE = table(NMSEvar,NMSEf,NMSElat,NMSEear,NMSElat2)

% Load additional data
% MOVE       
load movdat  %get Movm: column 1-dates, 2-values

dMovm = Movm(1+12:end,2) - Movm(1:end-12,2);
begm1fit1dMovm = begfit1<=Movm(:,1) & Movm(:,1)<=endfit1;
begm1fit2dMovm = begfit2<=Movm(:,1) & Movm(:,1)<=endfit1;

% CPI
load cpidat %get CPIr dCPIr
%make same lenght as rest 1979.10 to 2023.12
CPIr = CPIr(datenum(1979,10,1)<=dCPIr & dCPIr<=datenum(2023,12,31));
dCPIr = dCPIr(datenum(1979,10,1)<=dCPIr & dCPIr<=datenum(2023,12,31));
dCPI = log(CPIr(1+12:end)) - log(CPIr(1:end-12));

% Central Bank gold holdings IFS IMF: load only 1979.10 to 2023
CBR = xlsread('GoldResInternational_Financial_Statistics.xlsx', 1, 'BH185:VR185');
dCBR = log(CBR(1+12:end))' - log(CBR(1:end-12))';


% OLS regression of gold prices on other factors
% Marginal R-squared for each of the 6 variables for 4 sample perioeds
% Full sample
dUSDI = log(USDI(1:end-12)./USDI(13:end));


stats = regstats2NW12(dpdat, [dy10cal ] ,'linear'); 
            rsqad1=stats.adjrsquare; beta1f=stats.beta; pval1f=stats.hac.pval;
            
stats = regstats2NW12(dpmod, [dy10cal ] ,'linear'); 
            mrsqad1=stats.adjrsquare; mbeta1f=stats.beta; mpval1f=stats.hac.pval;
            
stats = regstats2NW12(dpdat, [dy10cal dUSDI ] ,'linear'); 
            rsqad2=stats.adjrsquare; 
stats = regstats2NW12(dpdat, [dy10cal dUSDI dy10cal-dy1cal] ,'linear'); 
            rsqad3=stats.adjrsquare; 
stats = regstats2NW12(dpdat, [dy10cal dUSDI dCPI ] ,'linear'); 
            rsqad5=stats.adjrsquare; 
stats = regstats2NW12(dpdat, [dy10cal dUSDI dCBR ] ,'linear'); 
            rsqad6=stats.adjrsquare; 
            
AddR2f = [rsqad1 rsqad2 rsqad3 nan rsqad5 rsqad6]';          

% Later 20 years
stats = regstats2NW12(dpdat(begm1fit1(13:end)), dy10cal(begm1fit1(13:end)) ,'linear'); 
            rsqad1=stats.adjrsquare; beta1lat=stats.beta; pval1lat=stats.hac.pval;
            
stats = regstats2NW12(dpmod(begm1fit1(13:end)), dy10cal(begm1fit1(13:end)) ,'linear'); 
            mrsqad1=stats.adjrsquare; mbeta1lat=stats.beta; mpval1lat=stats.hac.pval;            
            
stats = regstats2NW12(dpdat(begm1fit1(13:end)), [ dy10cal(begm1fit1(13:end)) dUSDI(begm1fit1(13:end))] ,'linear'); 
            rsqad2=stats.adjrsquare;            
stats = regstats2NW12(dpdat(begm1fit1(13:end)), [ dy10cal(begm1fit1(13:end)) dUSDI(begm1fit1(13:end)) dy10cal(begm1fit1(13:end))-dy1cal(begm1fit1(13:end))] ,'linear'); 
            rsqad3=stats.adjrsquare;                
stats = regstats2NW12(dpdat(begm1fit1(13:end)), [ dy10cal(begm1fit1(13:end)) dUSDI(begm1fit1(13:end)) dMovm(begm1fit1dMovm(13:end))] ,'linear'); 
            rsqad4=stats.adjrsquare;    
stats = regstats2NW12(dpdat(begm1fit1(13:end)), [ dy10cal(begm1fit1(13:end)) dUSDI(begm1fit1(13:end)) dCPI(begm1fit1(13:end))] ,'linear'); 
            rsqad5=stats.adjrsquare;  
stats = regstats2NW12(dpdat(begm1fit1(13:end)), [ dy10cal(begm1fit1(13:end)) dUSDI(begm1fit1(13:end)) dCBR(begm1fit1(13:end))] ,'linear'); 
            rsqad6=stats.adjrsquare;     
 
AddR2lat = [rsqad1 rsqad2 rsqad3 rsqad4 rsqad5 rsqad6]';            
            
 
% Early  20 years
stats = regstats2NW12(dpdat(begm1fit1(13:end)==0), dy10cal(begm1fit1(13:end)==0) ,'linear'); 
            rsqad1=stats.adjrsquare; beta1ear=stats.beta; pval1ear=stats.hac.pval;
            
stats = regstats2NW12(dpmod(begm1fit1(13:end)==0), dy10cal(begm1fit1(13:end)==0) ,'linear'); 
            mrsqad1=stats.adjrsquare; mbeta1ear=stats.beta; mpval1ear=stats.hac.pval;            
            
stats = regstats2NW12(dpdat(begm1fit1(13:end)==0), [ dy10cal(begm1fit1(13:end)==0) dUSDI(begm1fit1(13:end)==0)] ,'linear'); 
            rsqad2=stats.adjrsquare;            
stats = regstats2NW12(dpdat(begm1fit1(13:end)==0), [ dy10cal(begm1fit1(13:end)==0) dUSDI(begm1fit1(13:end)==0) dy10cal(begm1fit1(13:end)==0)-dy1cal(begm1fit1(13:end)==0)] ,'linear'); 
            rsqad3=stats.adjrsquare;                
 
stats = regstats2NW12(dpdat(begm1fit1(13:end)==0), [ dy10cal(begm1fit1(13:end)==0) dUSDI(begm1fit1(13:end)==0) dCPI(begm1fit1(13:end)==0)] ,'linear'); 
            rsqad5=stats.adjrsquare;  
            
stats = regstats2NW12(dpdat(begm1fit1(13:end)==0), [ dy10cal(begm1fit1(13:end)==0) dUSDI(begm1fit1(13:end)==0) dCBR(begm1fit1(13:end)==0)] ,'linear'); 
            rsqad6=stats.adjrsquare; 
            
AddR2ear = [rsqad1 rsqad2 rsqad3 nan  rsqad5 rsqad6]';   

% Later since 2007
stats = regstats2NW12(dpdat(begm1fit2(13:end)), dy10cal(begm1fit2(13:end)) ,'linear'); 
            rsqad1=stats.adjrsquare; beta1lat2=stats.beta; pval1lat2=stats.hac.pval;
            
stats = regstats2NW12(dpmod(begm1fit2(13:end)), dy10cal(begm1fit2(13:end)) ,'linear'); 
            mrsqad1=stats.adjrsquare; mbeta1lat2=stats.beta; mpval1lat2=stats.hac.pval;
            
stats = regstats2NW12(dpdat(begm1fit2(13:end)), [ dy10cal(begm1fit2(13:end)) dUSDI(begm1fit2(13:end))] ,'linear'); 
            rsqad2=stats.adjrsquare;            
stats = regstats2NW12(dpdat(begm1fit2(13:end)), [ dy10cal(begm1fit2(13:end)) dUSDI(begm1fit2(13:end)) dy10cal(begm1fit2(13:end))-dy1cal(begm1fit2(13:end))] ,'linear'); 
            rsqad3=stats.adjrsquare;                
stats = regstats2NW12(dpdat(begm1fit2(13:end)), [ dy10cal(begm1fit2(13:end)) dUSDI(begm1fit2(13:end)) dMovm(begm1fit2dMovm(13:end))] ,'linear'); 
            rsqad4=stats.adjrsquare;    
stats = regstats2NW12(dpdat(begm1fit2(13:end)), [ dy10cal(begm1fit2(13:end)) dUSDI(begm1fit2(13:end)) dCPI(begm1fit2(13:end))] ,'linear'); 
            rsqad5=stats.adjrsquare;
stats = regstats2NW12(dpdat(begm1fit2(13:end)), [ dy10cal(begm1fit2(13:end)) dUSDI(begm1fit2(13:end)) dCBR(begm1fit2(13:end))] ,'linear'); 
            rsqad6=stats.adjrsquare;             
 
AddR2lat2 = [rsqad1 rsqad2 rsqad3 rsqad4  rsqad5  rsqad6]';  

AddR2var = {'10-year real rate'; 'Y10&Real exchange rate'; 'Y10&Rer&Term spread 10-1'; 'Y10&Rer&MOVE';'Y10&Rer&CPI';'Y10&Rer&CBRes' };

Table6bottomAddR2 = table(AddR2var ,AddR2f ,AddR2lat ,AddR2ear ,AddR2lat2  )

          

[toc timesolveFxmodel]